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(p ! Abstract 

Establishing a direct link between individual based models and the corre- 
ct . spending population description is a common challenge in theoretical ecol- 
c/5 , ogy. Swarming is a prominent example, where collective effects arising from 

interactions of individuals are essential for the understanding of large-scale 

spatial population dynamics, and where both levels of modelling have been 

^ often employed without establishing this connection. 

r^ ■ Here, we consider a system of self-propelled agents with velocity alignment 

O ! in 2D and derive a mean-field theory from the microscopic dynamics via a 

nonlinear Fokker-Planck equation and a moment expansion of the probability 



density. We analyze the stationary solutions corresponding to macroscopic 
>^ ' collective motion (ordered state) and the disordered solution with no collec- 

^ ■ tive motion in the spatially homogeneous system. In particular, we discuss 

\Q . the impact of different propulsion functions governing individual dynamics. 

Our results predict a strong impact of individual dynamics on the mean field 

C^ I onset of collective motion (continuous vs discontinuous). In addition to the 

macroscopic density and velocity fields, we consider the effective "temper- 
ature" field, measuring velocity fluctuations around the mean velocity. We 
show that the temperature decreases strongly with increasing level of coUec- 

k> i tive motion despite constant fluctuations on individual level, which suggests 

H [ that extreme caution should be taken in deducing individual behavior, such 

as, state-dependent individual fluctuations from mean-field measurements 
[Yates et al. (2009) Proc. Natl. Acad. Sci. USA 106:5464-5469]. 
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1. Introduction 

Collective motion of living organisms, such as exhibited by bird flocks, fish 
schools or insect swarms is an ubiquitous and fascinating self- organization 
phenomenon, which has attracted scientists with very different backgrounds, 
ranging from biology and ecology to mathematics, physics and engineering. 

The understanding of grouping and collective behav ior is essential for 



our understanding of large scale dynamics of populations ( lOkubo and Levinl . 



200 ll : iKrause and Ruxtonl . |2002| ). The effective interaction of individuals 
leading to the onset of collective motion of up to thousands or millions of 
individuals, does not only have a strong impact on dispersal and migration 
of entire populations, but poses also interesting theoretical questions about 
modelling ar id the corresponding mathenaatical description of such ecologi- 



cal systems (jParrish and Edelstein-Keshetl . Il999l ) . The individual based - or 



"Langragian" - approach provides a natural framework for a detailed descrip- 
tion of individual dynamics in terms of (stochastic) equations of motion. The 
downside of the approach is the difficulty to obtain analytical results on the 
behavior of many interacting individuals. Thus, the approach often relies on 
extensive numerical simulations, which make it difficult to gain a profound 
understanding of the general dynamical behavior of the system. Further- 
more, despite the continuous progress in the development and application 
of novel computational techniques, large scale simulation of individual based 
model on ecological length and time scales are still extremely expensive in 
terms of computational time. Thus, at the level of populations, it seems 
more reasonable to take the "Eulerian" viewpoint, where the behavior of the 
ecological system is described by partial differential equations governing the 
dynamics of mean field observables such as the population density, which 
may even allow us to obtain analytical results. 

It is obvious that the different viewpoints are not independent: popula- 
tions consist of individuals; their dynamics origin from individual behavior. 
The question is, whether, for a given microscopic individual based model, 
a corresponding macroscopic or mean-field description can be derived. De- 
pending on the complexity of the individual based model this might be not 
feasible, but for simple models which focus on few essential mechanisms, it 
is possible to establish a direct link between the two levels of description. 



Here, starting from a simple individual based model of collective motion, 
we will derive systematically the corresponding mean-field equations. We 
consider a model of individuals interacting only via a velocity-alignment in- 
teraction, which tends to harmonize the speeds and directions of motion of 
neighboring individuals. The description of individuals dynamics is based 
on the concept of so-called Active Brownian particles, which can be de- 
scribed by stochastic equations of motion with a velocity- dependent friction 
function which, d e pending on the velocity, may also assume nega tive values 



(jSchweitzer et al.l . Il998l : lErdmann et al.l . l2000l : ISchweitzerl . 120031 ). Thus, in 



the following we will refer to individuals as "particles" or "agents" . 

Two types of velocity alignment can be distinguished: Particles with ne- 
matic interaction align either parallel or anti parallel, whereas a polar interac- 
tion acts strictly to wards parallel alignment of individual velocity vectors (see 



Ramaswamy 



2OIOI for a review). The velocity-alignment interaction studied 



here is polar (ICzirok et al.l . 119961) and can be seen as a continuous ve rsion of 
the well known Vicsek-model (IVicsek et al.l . Il995l : I Chat e et al.l . l2008l ). It re- 



duces for pairwise interaction of self-prope lled particles w i th con stant speed 
to the polar-alignment model studied by iPeruani et al.l ( 120081 ). Recently, 
there has been a number of publication on the coarse-grained description 
of self-propelled particles w ith velocity alignment. The ansatz proposed by 
Toner and Tul (Il998l . Il995l ) is based on the formulation of mesoscopic equa- 
tions of motion for the den sity and velo c ity fields using symmetry and conser- 



vation laws (see also, e.g ., iToner et al.l . l2005t iRamaswamyl . I2OIOI ). Recently, 



Bertin et al.l (120061 12009| ) derived hydrodynamic equations of interacting self- 
propelled particles by a Boltzmann approach. An alterna tive derivat ion of 
kinetic equations for the Vicsek model, was presented by llhld (I2OIII ). Ihle 
derives the corresponding transport coefficients by taking into account mul- 
tiparticle collisions and shows that the corresponding results are in good 
agreement with numerical simulations of the Vicsek model. Bertin et al. as 
well as Ihle assume in their derivations of coarse-grained description to be 
close to the critical point for the onset of collective motion (small center of 
mass velocity). 

Here, in contrast to the Vicsek model, we consider individual dynamics in 
terms of stochastic differential equation (Langevin equation) with continuous 
time. We do not assume a constant speed or restrict to the vicinity of the 
critical point for the onset of collective motion. The mean field equations are 
derived in a systematic way starting from the microscopic Langevin equations 
by formulating the corresponding Fokker-Planck equation for the probability 



density. Hereby, we consider directly multi-individual interactions as each 
individual couples to the mean velocity of its neighbours. In addition to the 
density and velocity fields, we consider explicitly the effective temperature 
field of the active Brownian particle gas. In this work, we focus on the impact 
of different (nonlinear) velocity-dependent friction functions on the onset of 
collective motion. 

Although the kinetic theory derived in this paper, can in principle be 
used to analyze the spatially inhomogeneous solutions, we restrict ourselves 
here to the discussion of the spatially homogeneous case. Please note that 
recent publications show the instability of the homogeneous solut i on in large 
syste r ns w i th sho rt-ranged velocity alignment (JBertin et al.l . l2006l l2009l : iLed . 
20 id : llhld . I2OIII ). Thus, the results presented here can not be generalized to 
the so-called thermodynamic limit with local interactions. 

The derivation of the kinetic theory is based on the formulation of mo- 
ment equations of the corresponding probability distribution. This approach 



has been employed, for example, by iRiethmiiller et al 



(119971) to analyze the 



behavior of a quasi one-dimensional granular system. lErdmannI (J2003[ ) used 
it to analyze the mean field of non -interacting active particles with Rayleigh- 
Helmholtz friction. Only recently, iRomanczuk and ErdmannI (120101 ) used the 
approach to analyze the collective motion of active particles in one spatial 
dimension. 

In general, for a system far from equilibrium the probability distribution is 
not Gaussiar i and a correct descri ption requires infinitely many moments (see 



for example iPawulal . Il967l . Il987l ). Thus, depending on the detailed model. 



it may be necessary to neglect higher moments in order to obtain a closure 
of the system of moment equations. The approximation of a non-Gaussian 
probability distribution by a finite number of moments may lead to unphys- 
ical behavior, such as negative values or artificial oscillations of the (approx- 
imated) probability distribution. Therefore, we will compare our analytical 
results to numerical simulations of the microscopic system. 



2. Active Brownian particles with velocity alignment 

We consider a system of A^ individuals with mass m in two spatial dimen- 
sions. For simplicity, we neglect the finite size of individuals and consider 
them to be point-like active Brownian particles. The evolution of the particle 
positions Tj = rj(t) and velocities Vj = Vj(t) is described by the following set 



y 




X 



Figure 1: Scheme of the velocity-ahgnment interaction. The focal individual (blue/dark 
grey) interacts with all neighbors (red/bright grey) within the metric interaction range e. 



of stochastic equations of motion (i = 1 . . . A^): 



(1) 
(2) 



The first term on the right hand side of Eq. ([2]) is a friction/propulsion 
force, which describes the deterministic velocity dynamics of non-interacting 
active particles. This velocity-dependent friction is negative at low speeds, 
which leads to an acceleration of the individual. As a consequence an indi- 
vidual has a preferred speed different from zero. In the following sections, 
we will discuss active particles with two different velocity -dependent fric- 
tion functions: the n onlinear Rayleigh-Helmholtz friction ( jRayleighl . Il894 : 
Erdmann et al.l. 120001). as well as a variant of the linear Sc hienbein-Gruler 



friction (jSchienbein and Grulerl . Il993l : lErdmann et al.l . 12000 ). 



T he second term is a velocity-alignrn ent interaction (JNiwal . Il994j : lOkubo and Levin 



200ll : iRomanczuk and Erdmanru . |2010| ) , with 



N,„ 



U£,i = Ue,i(ri,t) 



iV. 






(3) 



being the mean velocity of particles within the neighborhood Se,i defined via 



the radius e > around the focal particlcl: Fj G S'^^j if |rj — r^l < e. The 
ahgnment acts towards harmonization of the velocity of the focal particle 
with the local average velocity of its neighbors. The alignment strength 
fj, = 1/Ta > determines the relaxation time Ta of the velocity of the focal 
particle towards the average velocity of surrounding particles (see Fig. [1]). 
For solitary particles, with no neighbors, or a system in a perfectly ordered 
state where all particles move with exactly the same velocity, the alignment 
force vanishes as u^ j = Vj. On the other hand, for a large number of neighbors 
Nfr^i ^ 1 moving with random velocities (disordered state), the mean velocity 
vanishes u^^j = and the velocity alignment force leads to an additional 
"social" friction — /iVj. 

The last term on the right hand side of Eq. ([2]) accounts for is stochastic 
force, which introduces a random component into the motion of each individ- 
ual. For each individual it is given by a Gaussian random force with intensity 
D independent on the current velocity of the focal individual as well as on 
the motion of other individuals. Here, ^(t) is a Cartesian vector with un- 
correlated components given by a stochastic Gaussian process: {C,k{t)) = 0, 

m)ut')) = SMS{t-t'). 

Throughout this work we will consider the motion of individuals in a 
rectangular spatial domain of the size L^ with periodic boundary condition 
(torus). Furthermore we rescale all forces by the mass of individuals, by 
setting m = 1. 

We introduce the probability distribution P(r, v, t) dr dv, which deter- 
mines the probability to find an individual at time t, in the spatial region 
[r, r -|- dr] moving with velocity within the interval [v, v + civ] . In the fol- 
lowing, we make the mean field assumption for the alignment term N^^i ^ 1 
and derive the mean field theory corresponding to the microscopic dynam- 
ics given in Eqs. ([T]), ([2]) via the formulation of moment equations of the 
corresponding probability distribution. We define the n-th moment of the 
fc-component of the velocity vector Vk = ve/c, with e^ being a canonical basis 
unit vector as 

(Vk) (r, t) = -^ I d^vlPiv, V, t) , (4) 

with k = x,y and n = 0, 1, 2, . . . . The spatial density of individuals p(r, v, t) 



^Please note that the sum includes the velocity of the focal particle v^. Thus, for a 
solitary particle N^ = 1 and u^^ = v^. 
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corresponds to the zeroth moment (n = 0) and gives us the normahzation of 
the probabihty density with respect to integration over the velocity space: 

p(r,t) = y"t/vP(r,v,t). (5) 

Multiplying the n-th moment by the density and taking the derivative with 
respect to time, we obtain the dynamics of the velocity moments: 

|(pK"» = /<iv„jf. (6) 

For more than one dimension the above definition can easily be extended to 
mixed moments as, for example, the covariance: 

(«) = - I dwvy^P{v,w,t) n,m > 0. (7) 

The Fokker-Planck equation, which determines the evolution of the dis- 
tribution function of individuals P, for the microscopic dynamics given in 
Eq. (ED reads 

^P(r, V, t) = -vVrP - Vv { [-7(v) V + ^ (u, - v)] P} + D^P . (8) 

The local average velocity u^ = Ue(r, v,t) in the alignment term depends 
implicitly on the probability density P, thus the above equ ation is a nonlinear 



Fokker-Planck equation as discussed, e.g., by lFrankl (120051 ). In the continuous 



description we may express u^ as an integral over the distribution function: 
u.(r, t) = ^ I dv' I dVVPir', V, t) . (9) 

Se represents the spatial neighbourhood of the position r, defined via a metric 
distance: r' G S'e if |r — r'| < e. 

Please note that in the limit £ — )■ 0, the local average velocity in the align- 
ment term u^ reduces directly to the mean field velocity u(r, t) = {{vx), (vy)), 
whereas for a finite value of e it can be seen as an approximation of u under 
the assumption of a constant density on the length-scale corresponding to e. 



3. Rayleigh-Helmholtz friction function 



At first we derive the mean field tlieory for active Brownian particles with 
the so called Rayleigh-Helmholtz friction, which reads 



-7(v) 



(a - /3v2)v, 



(10) 



with a,/3 > 0. The above friction (or propulsion) function in the has been 



studi e d in detail in the context of ac t ive Br ownian motion (lErdmann et al. 
2 pool . I2OO2I : ISchimanskv-Geier et al.l. 120051) . It was used in general mod- 



els of swarm dynamics (IMikhailov and Zanettd. Il999l: iRappel et al.l. Il999l: 



Erdmann et al.l . l2005l : iD'Orsogna et al.l . 120061 : lEbeling and Schimanskv-Geier 



20081 ) as well as for modeling of fi sh schools (JNiwal . Il994j. 119961). bacteria l 
swarming (JRomanczuk et al.l . 120081 ) and locusts swarms (JBazazi et al.l . I2OIOI ). 
The term av in f lTOj) leads to a acceleration of individuals at low velocities 
in the direction of motion, whereas the term — /3v^v represents a nonlinear- 
friction. The deterministic, stationary speed of individuals can be calculated 
from the balance of the acceleration and fr iction to Vn = \/a/l3. 



In analogy to the one- dimensional case (JRomanczuk and Erdmannl . I2OIOI ) 
we insert Eq. (|8]) with Eq. ( ITOj) into Eq. ([6]) and obtain for the n-th moment 

of Vr 



|(.«» 



V^ \ —Vxi;, uy 



d d 

dy 



_d_ 
d 



[- {a - P{vl + V^)) Vx + ^l{Ue,x - V:^)] 



-^['i^- /^(^^ + ^y)) ^y + ^("^>2/ ~ ^y)] 



+ D 









>Pdv. 



(11) 

Since the probability distribution approaches zero at infinity limy^^±oo P = 
0, the terms with derivatives with respect to Vy vanish. The terms with 
derivatives with respect to v^ can be partially integrated and we obtain 

|(pK))--|(p«-))-|(pK^.)) 



+ n {n-l) Dp{v:-'). 



(12) 



Please note that the dynamics of the n-th moment depend on the n+ 1-th 
moment. Thus, we obtain an infinite hierarchy of coupled moment equations, 
for which we have to define a closure condition. 

We rewrite the velocity vector as a sum: v = u + 6v, with a mean-field 
velocity u = {ux,Uy) and a vector of symmetric deviations around the mean 
5v = {6vx, Svy) with (5f^) = for odd a. . Furthermore, we assume that the 
deviations in x and y are independent, thus 

{Sv:Sv^) = {Sv:){5v^). (13) 

We have verified this assumption by numerical simulations, which show neg- 
ligible mean-field cross-correlations of velocity deviations. This can be also 
justified by symmetry considerations of the two possible stationary states 
discussed in detail further below: the completely isotropic disordered state 
with vanishing mean velocity |u| = 0, and the ordered state with |u| > 0, 
where the cross-correlations of the deviations parallel and perpendicular to 
the mean field velocity must vanish. 

From the above assumption we obtain a diagonal covariance matrix: 

^- - ^ {5vy5v,x) {Svp )-[o tJ- ^^^' 

In analogy to kinetic gas theory we will refer to the "vector" of the non- 
vanishing diagonal elements T = (Tj.,Ty) as temperature. Please note that, 
in contrast to simple gas at thermal equilibrium the components of T may 
differ at non-equilibrium steady state: T^ ^ Ty. 

Using u and T we may write the moments {v"^) as: 

(ffc) = Uk (15a) 

{vl)=ul + n (15b) 

(^f ) =ul + ?>UkTx (15c) 

{vl)=ul + Quln + Tl + ek^ (15d) 

with k = x,y and 6 u being the mean field temperature fl uctuations in k- 
direction (Erdmannj, l2003l : iRomanczuk and Erdmannl . 120101 ) defined as: 



Ok = {{{vk - Ukf - nf) = {6vt) - Tl (16) 

These higher order fiuctuations cannot be expressed by the mean velocity 
and temperature in the case of non-Gaussian distributions. The equations 
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governing the dynamics oi = {6x,dy) can be derived for nonlinear friction 
functions by considering the evolution of higher moments {n > 2). 

For the mixed moments {v^v^") we obtain the following expressions: 

{Vx Vy) = Ux Uy (17a) 

{vx fy) =Ux ul + Ux Ty (17b) 

{Vl Vy) = Ul Uy + Uy Tx (l7c) 

(^x ^l) =ulul + ul Ty + ul Tx + Tx Ty. (17d) 

We insert Eqs. ( !T5|) and (TT7|) in the equation for the moment dynamics ( !T2|) 
and, by carrying out the calculations up to the second order (n = 2), we 
arrive at the following set of partial differential equations 

-p = - Vr (pu) , (18a) 

^ + uVrMx = <yUx - (3Ux (u^ + 3Tj, + Ty) + li{Ue,x - U^) 



dt 

dx p dx' 

1 ron 



(18b) 



2\ dt 



+ uVrTj ={a- i2)Tx - l3Tx (u^ + 2ul + T, + Ty) 



^D-T.^. (180) 

Here, we have only given the equation for x-components of the mean velocity 
and temperature; the corresponding equations for the y-component can be 
directly obtained by interchanging x and y. 

In the following we will focus on the spatially homogeneous system, which 
corresponds to the case of global coupling. In this case the above equations 
simplify to ordinary differential equations for the mean field velocity and 
temperature (x-component): 

— ^ =aux - (5ux (u^ + STj. + Ty) , (19a) 

au„ - (5uy [v? + Tx + "iTy) , (19b) 



dt 



Uy 



1 y^T' 

" {a-p)Tx-/3Tx{u^ + 2ul + Tx + Ty)-f3ex + D, (19c) 



2 dt 



1 J'T"' 

' -{a - p)Ty - f3Ty (u^ + 2ul + T, + Ty) - (iOy + D. (19d) 



2 dt 
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We perform a closure of the moment equations by setting 6x,y = and 
obtain a self-consistent set of ODE's. This is a reasonable approximation at 
low noise intensities D. In general, finite higher order fluctuations dx,y > 
have an inhibiting effect on the temperature T, which for sufficiently large 
D may be not negligible. 

We can further simplify the above set of equations by choosing a reference 
frame where u^ = u\\ = u corresponds to the mean field velocity, whereas the 
orthogonal component vanishes Uy = u± = 0: 



du 

Hi, 
IrfTii 



-au 



2 dt 
IdTi 



2 dt 



(3u {u^ + 3T|| + Ti_) (20a) 

a - ii)T\\ - I3T\\ (3^2 + T|| +Ti_)+D (20b) 

a - ii)T^ - (3Ti_ (m^ + Ty +Ta_)+D (20c) 



with T|| and T± being the temperature components parallel and perpendicular 
to the mean field direction of motion. In the stationary disordered state 
(m = 0), the temperature components can be easily calculated from fl20|) and 
the corresponding solution reads 

Ml =0, (21a) 

lli = l±,i = li= j^ . (2ibj 

In the case of vanishing noise D = 0, the ordered solution can be immediately 
obtained as u = \fotJ^ and T\\=T^ = 0. For D > 0, it is evident that the 
temperature component parallel to the direction of motion is smaller than 
the perpendicular one: T\\ < T±. 

For the general ordered state, with u > and -D > 0, we were so far not 
able to obtain explicit stationary solution for u, Ty and T± of the above ODE 
system (1201) but the stable and unstable solutions can be determined by a 
numerical continuation metho ds as, for exaraple, provided by t he numerical 



software XPPAUT/AUT097 flDoedei Il98ll : lErmentroutl . 120021) 



A possible ansatz to find an explicit solution is the reduction of the di- 
mensions in the problem: We use the fact that at a fixed time t we may 
always find a coordinate frame where Ux = Uy = u. In this coordinate frame 
due to the symmetry of the involved equations we obtain also T^ = Ty = T. 
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Figure 2: Comparison of the mean field speed |m| obtained from Langevin simulations 
(symbols) of the RH-model in two spatial dimensions with the results of the mean field 
theory for the homogeneous case. The black lines represent the solutions obtained from 
the full system of mean field ODE's. The red lines represent the mean field solutions from 
the reduced system. The stable solutions are shown as solid lines, whereas dashed lines 
indicate the unstable solution. The simulations were performed with periodic boundary 
condition and with the disordered state as initial condition. Other parameters: a = 1, 
I3 = I^L = 200, e = 20, iV = 4096. 
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noise D 



0.01 0.1 
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Figure 3: Numerical results for stationary temperature components, Tn (blue/dark gray 
circles) and Tj_ (black squares), obtained from Langevin simulations of the RH-model with 
mean field predictions for weak velocity alignment (left) and strong velocity alignment 
(right). The solid (dashed) lines represent stable (unstable) solutions. The solutions of 
the full system for Tn {T±) are indicated by blue/dark gray (black) lines. The red lines 
show the mean field solutions from the reduced system. The simulations were performed 
with periodic boundary condition and with the disordered state as initial condition. Other 
parameters: a = 1, /3 = 1, L = 200, e = 20, iV = 4096. 
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With this ansatz we reduce the original system f ll9p from a four dimensional 
system of ODE's to a two dimensional system in u and T: 



^=au-Pui2u^ + AT 



IdT 
2ltt 



{a-fi)T~l3T{Au' + 2T 



D 



(22a) 
(22b) 



This gives us a system of equations similar to the one-dimensional case, where 
we can analytically determine the stationary solutions for du/dt = dT/dt = 
to: 



«1,2 = 0, 



Ti, 



U3A 



n = T^ 



^5,6 



n = n 



a — fi± 


v/(a - /i)2 + 8/3D 


4/? 


J2a- 


-fi+y/{a + 


/i)2 - 24/3D 




VW 




a + fi - 


y/{a + fiy- 


-24/3D 


12/3 


J2a- 


- fi- v^(a + 


/i)2 - 24/3D 




VOP 




a + /i + 


Via+f^r- 


-24/3D 



12/3 



(23a) 
(23b) 

(23c) 
(23d) 

(23e) 
(23f) 



The type of the different solutions is completely analogous to the one-dimensional 
system: 



Ml,2, Ti 

"^3,4; ^3,4 
■^5,6; ^5,6 



disordered solution 
stable ordered solutions 
unstable ordered solutions 



The structure of the bif urcation diagram of the reduced sy stem fj22|) for d = 2 
is the same as for d = 1 ( Romanczuk and Erdmannl . 120101 ) with different crit- 



ical noise intensities determining the stability of the disordered and ordered 
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solutions: 

D 



(2d) _ a(2/i-a) 

d,cr\t o /P ' 



^/3 

(rv -I- 

D 



(2d) _ (Ct + Z^)^ 
'o,crit ~ 



24/3 

The three different regimes obtained from the reduced system can be sum- 
marized as: 

1. weak ahgnment (0 < ;U < ^): 

(a) bistabihty for weak fluctuations D < Do^cnt 

(b) only disordered solution stable for strong fluctuations D > -Do,crit 

2. intermediate alignment (| < /i < 2a): 

(a) only ordered solution stable for weak fluctuations D < -D^^crit 

(b) bistabihty for intermediate noise -Dc^ crit < D < Z^ccrit 

(c) only disordered solution stable for strong fluctuations D > -Do.crit 

3. strong alignment (/i > 2a): 

(a) only ordered solution stable for weak fluctuations D < Z^o.crit 

(b) only disordered solution stable for strong fluctuations D > D^crit 

The comparison of the stationary solutions of the reduced systems with the 
corresponding solutions of the full systems obtained with XPPAUT/AUTO 
reveals differences at low velocity-alignment strengths /x (Fig. I2]l3l) . The ve- 
locity of the stable ordered solution of the full system decreases more strongly 
and exhibits an earlier breakdown with increasing D. Furthermore, from the 
position of the disordered branch, it can be deduced that the basin of at- 
traction of the ordered state for the full system at low noise D is larger than 
for the reduced two dimensional system. At large /i the differences between 
the two types of mean field solution vanish and the reduced system gives a 
good approximation as shown in Fig. |2l The reason for the discrepancy be- 
tween the two mean field solutions at low /i is the fact that the reduction of 
the system dimensions "throws away" all information about the asymmetry 
of temperature components parallel and perpendicular to the mean velocity. 
But at large /x both temperature coefficients are dominated by the — /xT^ term, 
so that the asymmetry in the temperature components becomes negligible 
as shown in Fig. [31 Please note, that in the two-dimensional phase-space 
projections used in Fig. [3l distinct non-overlapping solution with \u\ = 
and |u| > appear on top of each other. 
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In general, without knowing the temperatures Ty and Tj_, the mean speed 
as the order parameter can be written as 



\u\ = ^vl - 3Tj| - T^ (24) 

here we used Vq = a/ (3. In the limit of large // close to the critical noise, where 
a, (3 <t^ fi, D, we may approximate the temperature as Ty = T± = T = D/fi 
and we obtain a simple expression for the ordered state 




u\ = \vi . (25) 



In this limit, the critical noise can be approximated as -Dd,crit ~ "^0/^/4 = -Dcrit 
and the above equation may be rewritten as: 

|M|=2/i-5(D,,i,-D)l, (26) 

which is the standard form of the order parameter for a continuous (second 
order) phase transition. 

In order to check the analytical results, we have performed large scale 
Langevin simulation of the microscopic system with periodic boundary con- 
dition. In two dimensions, the direction of motion in the ordered state u/|u| 
can freely diffuse. Due to always present fluctuations in a finite system the 
direction of motion changes in time. Thus, the ordered state is characterized 
in simulations through the non-vanishing mean speed of the system: 



(|u|) 



A 




(27) 



Here, N is the total number of individuals and (■) indicates the temporal 
average after the system reached the stationary state. 

The results of microscopic Langevin simulations in two spatial dimensions 
d = 2 aX high densities and L/e = 10 are in a good agreement with the 
semi-analytical results for the full mean field system at low n and with both 
solutions types at large fi. The temperature components Tj| and T± differ 
significantly at low fi and the solution of the reduced system is not able to 
describe the system correctly. But for large /i the differences between the 
different components become negligible and the different stationary mean 
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field solution and the numerical results collapse practically on a single line 

(Figs. El ED- 

Even if the corresponding basin of attraction of the disordered state is re- 
duced in comparison to the one-dimensional case (JRomanczuk and Erdmannl . 
20101), the mean field theory predicts nevertheless its stability at low /x and 
low noise D. But in contrast to the one-dimensional case, we were not able 
to observe a stable stationary disordered state in numerical simulations even 
at very low /i and low but non- vanishing D (Fig. [2]). This is a consequence 
of neglecting all higher order fluctuations 9k in our mean field ansatz. In 
two dimensions it has a stronger impact on the stability of the disordered 
state than in the one-dimensional case in the respective parameter region. 
In two spatial dimensions at low /i and low D there is no energetic barrier 
between different directions of motion. Here, the individuals can change their 
direction of motion by continuous angular drift or diffusion. Thus, any finite 
fluctuation in u at vanishing noise will be amplified and eveii t ually will lead 
to perfect alignment (see also lEbeling and Schimansky-Geierl . |2008| ). 

Recently, it was shown that the spatially hom ogeneous state is unstable 



in systems of interactin g self-pr op elled parti cles (iBertin et al.l . |2006| . 12009 



Simha and Ramaswamyl . l2002bl Jal: llhld . I2OIII ). Thus, the assumption of a 
spatially homogeneity is only an approximation, which is justified only for 
sufficiently large e (upper panel. Fig. H]). For L ^ e strong density inhomo- 
geneities appear, such as travelling bands (lower panel. Fig. Hj), which affect 
the global behavior of the system. 



4. Schienbein-Gruler friction function 

As a second example, we consider a system of active Brownian p articles 
with linear speed dynamics studied by ISchienbein and Grulerl (19931). Here 



we us e a corresponding friction function as introduced by lErdmann et al 
(j2000l ) in two spatial dimensions: 



-7(^vjv 



-a 



1-^ 



(28) 



where, Vq is the preferred speed of individuals (stationary speed), and a 
is now the inverse relaxation time of the velocity dynamics. Please note 
that the above friction/propulsion term, shows a discontinuity at the origin 
V = 0, which may be eliminated by taking into account the possibility of 
backwards motion of individuals with respect to their heading (JRomanczukl . 
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Figure 4: Snapshots of the ordered state of active Brownian particles with RH-friction and 
velocity alignment for two different values of e = 1.0 (left) and 5.0 (right). The large arrow 
indicates the mean direction of motion. (Other parameter values: N — 8192, /i = 1.0, 
D = 0.1, L = 100) 



20 111 ) neglected by ISchienbein and Grulerl (119931 ) . But the much most im- 
portant property of the above friction function, irrespective of the details at 
vanishing speed, is its linear behavior around the stationary speed Vq. 

Following the same approach as in the previous examples we arrive at the 
following moment equation for the x-component: 



l'^<' 



- ^ (p {<*')) -^(p {< %}) 



+ n p 



iC 



a 



v^^ 



K) +M«.,x.«-^)-K)) 



+ n (n-l) Dp{v:-') 



(29) 



The above equation differ from the one obtained in the Rayleigh-Helmholtz 
case ( TT2I) only in terms associated with the friction function. Due to the 
absolute value of the velocity |v| in ( l28l) . we obtain in the moment evo- 
lution equation not only linear combination of moments but also the term 
vl + f^), which we have to take care of. Here, we make the following 
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approximation : 




(v/^^F+^> ^ui + ul + T, + Ty 



(30) 



This approximation, results in a closure of the final mean field equations. 
For clarity, we will use the following abbreviation for the denominator in the 
above expression 



Vt := VAvL, T) = Jul + ul + T,. + T,. 



Thus, we obtain the following system of mean field equations for the /cth- 
component of the mean field velocity and temperature {k = x,y): 

duk ^ f (y.VQ \ , dTk Tk dp 

+ uVrUfc = — a]uk + p[u^,k - Uk) - -^ ^-, (31a) 



dt \^T ) ' dx p dx' 

+ uVrTfe = — a-p]Tk + D- Tk^r-, (31b) 



which together with the continuity equation (I18ap determine the evolution 
of the mean field system. 

By considering again the simplest case of a spatially homogeneous system, 
we arrive at the following set of ordinary differential equations [k = x,y or 

duk favo \ 

- a ] Uk, (32a) 



dt \Vt 

1 dTk f OiVQ 



a-p]Tk + D. (32b) 



2 dt \Vq 

From the symmetry of the temperature equations it follows directly Tj. = Ty 
or by choosing the appropriate reference frame Ty = T±. Using the mean ve- 
locity squared u^ = u'^+Uy as a. variable instead of the individual components 
{ux, Uy) the above equations simplify to: 



—u^ =2a {^-l]u\ (33a) 

-^ = ^-a-Hr, + D. (33b) 
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There are two stationary solution of Eq. 
components reads: 



lal in terms of the temperature 



u 



u 



0, 

^0 - T\ 



Ti 



(34a) 
(34b) 



The full system of equations has in total three stationary solution: a single 
ordered state solution with \u\ = yv? > and two disordered solutions with 
0. The single ordered solution reads: 



\u 



w 1 



2D 



T -^ 



whereas two stationary disordered solutions are given as: 

|m|2,3 = 0, 



Ti 



D 



k,2 



+ 



a + /i \2(a + /i) 



T> 



D 



fc,3 



+ 



avo 



a + n \2{a + jj,) 




8D{a + fi) 



2„2 



a^v. 



8D{a + fj.) 



2,, 2 



a^v, 







(35) 

(36) 
(37) 

(38) 



We identify Tk,2 with the temperature of the disordered phase. The second 
disordered solution Tk^s is not considered, as it is always less stable than Tk,2 
and yields unphysical result of vanishing temperature of the disordered phase 
in the limit -D,/i — )■ 0. For noninteracting self-propelled particles, for which 
only the disordered state exists, the temperature has to coincide with the 
average kinetic energy per particle: Vq/2. It can be easily shown that Tk^2 
satisfies this condition in the respective limit. 

For noise intensities below the critical noise intensity 



D, 



crit 



^oV/2, 



(39) 



the homogeneous system is in the ordered state {\u\ = \u\i, T^ = Tk^i), 
whereas above the critical noise intensity no collective motion takes place 
(|m| = \u\2, Tfc = Tfc2)- For vanishing noise, D — ?■ 0, the temperature com- 
ponents Tfc vanish and the mean field speed is given as |m| = vq. For t>o = 0, 
the ordered state is always unstable and the temperature of the disordered 



19 



CD 
CD 
ft 

a; 



1.25 



0.75 



0.5 



0.25 









. -*<u^ 


tt-H. 


- 


\* 


\ \ 


- 


. ^ = 0.05 \ 

■ tJ. = 5.0 

1 


\ \ 


- 



2.5 



0.001 0.01 



0.1 



10 






0.5 




noise D 



0.001 0.01 0.1 
noise D 



Figure 5: Comparison of the simulation and analytical results for the mean field speed |u| 
(left) and temperature |T| (right) versus noise strength D for the Schienbein-Gruler fric- 
tion. The results of Langevin simulations are shown as symbols. The solid lines represent 
the stationary solutions according to Eq. (j41ap . The simulations were performed with 
periodic boundary condition and with the disordered state as initial condition. Simulation 
parameters: a = 1.0,wo = 1-0, £ = 10.0, L ^ 100.0, N ^ 16384. 
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Figure 6: Comparison of the simulation and analytical results for the mean field speed |u| 
(left) and mean temperature |T| (right) versus velocity alignment /i for the Schienbein- 
Gruler friction. The results of Langevin simulations are shown as symbols. The solid 
lines represent the stationary solutions according to Eq. ()41ap . The simulations were per- 



formed with periodic boundary condition and with the disordered state as initial condition. 
Simulation parameters: a = 1.0, uo = 1-0, e = 10.0, L = 100.0, N = 16384. 
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Figure 7: Left: Coniparison between |m| from Eq. P5)) . |u| from Eq. ([55]) and nmnerical 
results for the Schienbein-Gruler friction; Right: The corresponding temperature compo- 
nents T|| = D/{a + ^) (blue/dark gray line), T± = D/ ji (red/light gray line) and Tk^2 from 
([37l) (black line) together with results of numerical simulations (symbols). The Langevin 
simulations were performed with periodic boundary condition and global coupling. Simu- 
lation parameters: a = 10.0,uo — 1.0, N — 16384. 



state reduces to T^ = D /[a + /i). Inserting the critical noise intensity into 
( 135|) yields for the mean speed 



\u\ 



2^fi-^(D-Dc) 



(40) 



This corresponds to a second-order phase transition, not only in the limit of 
large /i, as in the case of Rayleigh-Helmholtz friction. 

The full stationary solution of the homogeneous system may be written 
as: 



\u\ 



n 



\u\i for D < Dcrit 

for D>D,,,t 

Tk,i for D < Dcrit 

Tfc ,2 for D > Dcrit 



(41a) 
(41b) 



As Tfc^i = Tfc_2 = "^0/2 for D = Dcrit the temperature is continuous but in gen- 
eral not continuously differentiable with respect to the bifurcation parameter 
(e.g. D, fi). 

In order to test the theoretical predictions, we have performed large scale 
Langevin simulation of the microscopic system either with global coupling 
6 > L, or with local coupling e = 1/lOL. At high densities, p = {Ne^)/ L"^ 3> 
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1, both situations yield similar results due to the large correlations length in 
the ordered phase with /corr ^ ^- Thus, also for local coupling the system 
can be assumed as spatially homogeneous, if the interaction range e is not 
to small with respect to the system size L. All simulations were performed 
with periodic boundary conditions and with random initial positions and 
velocities (disordered state). 

As shown in Figs. and the numerical results for the mean speed and 
the total temperature 



are in good agreement with the theoretical prediction. But a closer look re- 
veals some systematic deviations. The theoretical results for the temperature 
T are larger than the ones obtained in the simulations in particular at low fi. 
The effect appears for both ordered and disordered state as well as for global 
and local coupling. This systematic difference vanishes for uq = 0, thus it 
may be associated with the moment approximation in the corresponding term 
of the friction function (Eq. fl30|) ). Furthermore, as in the case of Rayleigh- 
Helmholtz friction, the temperature components parallel and perpendicular 
to the direction of motion are not equal, in contrast to the predictions of the 
mean field theory. Whereas the perpendicular component agrees well with 
the theoretical prediction of T±_ = D/n, we observe in general smaller values 
of the parallel component. 

Starting from the perfectly ordered state in the homogeneous state, it can 
be shown that for small noise intensities the components of the temperature 
can be approximated as: 

- « f,^£. (42) 



a + /i /i 

Using Eq. I34bt we obtain then an alternative solution for the mean field 
speed: 



In the limit of small a or large fi, this result converges towards the original 
mean field result. This solution {\u\, Ty, T±) for the ordered state yields a 
better agreement with numerical results not to close to the critical point but 
deviates from the numerical solutions close to the onset of collective motion 
(see Fig. [7]). 
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5. Conclusion and Outlook 

We have shown how starting from the microscopic dynamics of active 
Brownian particles with velocity alignment a corresponding mean field theory 
for the coarse grained density p(r, t), velocity u(r, t) and temperature T(r, t) 
can systematically be derived. 

In general, depending on the details of the individual dynamics, approx- 
imations have to be applied in order to obtain a self-consistent set of equa- 
tions. Here, we have discussed two different friction/propulsion functions 
governing the motion of individuals. We have compared the results of the 
mean field theory obtained under the assumption of a spatially homogeneous 
case, and have shown that they are in good agreement with numerics for a 
wide range of parameters. On the other hand, the approximations made in 
the derivation of the mean field theory, may result for some cases in deviations 
between the analytical predictions and numerics, which can be understood 
from a detailed analysis. 

The mean field approach provides important insights into collective dy- 
namics due to velocity alignment. In the case of a spatially homogeneous 
system (e.g. global coupling), where the local velocity field the individu- 
als are sensing, equals to the mean field velocity u^ = u the alignment force 
term vanishes in the velocity equation. Thus, in this limiting case the velocity 
alignment acts only on the effective temperature of the active Brownian par- 
ticle gas. It suppresses the individuals fluctuations around the current mean 
velocity leading to the stabilization of the ordered state with finite mean fi eld 



velocity (lEbeling and Schimanskv-Geieii . l2008l : lRomanczuk and Erdmanru . l2010l ) 



Recently, lYates et al.l (120091 ) analyzed the mean velocity and its fluctua- 
tion in the collective motion of locusts. Based on the observed dependence of 
mean velocity fluctuations they suggested that this implies state dependent 
fluctuations (multiplicative noise) on the level of individuals. In our case the 
"temperature" |T| is a measure of fluctuations of individual velocities around 
the mean. Although it does not correspond directly to the fluctuation ob- 
servable studied by Yates et al. in a finite system it is nevertheless related. 
The stationary values show a strong dependence of |T| on the correspond- 
ing mean velocity - increasing |T| with decreasing! u| - even in the absence 
of any multiplicative noise in the dynamics of individuals. This is a conse- 
quence of the nonlinear Fokker-Planck equation governing the dynamics of 
the probability distribution and the resulting effective state-dependent ve- 
locity potential. Thus our results suggest that in biological systems extreme 
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caution should be taken when deducing individual behavior of interacting 
agents based on mean field measurements, and that there might exist other 
alternative explanation of the o bserved mean field dynamics of locusts than 



the suggested "inherent noise" (JYates et al.l . 120091 ) . 



Furthermore, we have shown that for the homogeneous case in the limit of 
large coupling fi, where the impact of the individual speed dynamics (friction 
function) in the temperature equations becomes negligible, the scaling of the 
mean speed as order parameter close to the critical noise reads: 

\u\-[-^^) . (44) 

Thus, in this limit the onset of collective motion in a spatially homogeneous 
system takes place via a continuous (second-order) phase transition, irre- 
spective on the details of the velocity dynamics. This result agrees also with 
previo us results obtained fo r a system of self-propelled particles with constant 



speed (jPeruani et al.l . l2008l ) 



In one dimension it was shown, that for weak coupling (low fi) and a non- 
linear Reyleigh-Helmhol tz friction function, bistability of the ordered and dis- 



order ed state is possible (JMikhailov and Zanettel . ll999l : lRomanczuk and Erdmann 



2010l ). The corresponding mean field theory predicts a discontinuous tran- 
sition from an initially ordered state of collective motion to the disordered 
state with increasing noise intensity for nonlinear friction functions. But 
the results of direct large-scale numerical simulations of the individual based 
model reveal deviations from the predicted behavior. 

A possible explanation could be an extremely small basin of attraction 
of the disordered solution in the respective parameter region and the always 
present finite fiuctuations in Langevin simulations (largest system studied 
with global coupling: A^ = 32768). But a more probable cause, is the ap- 
proximation of the non-Gaussian probability distribution by a finite num- 
ber of moments, and the corresponding unresolved impact of higher order 
fiuctuations on the stability and bifurcation behavior of the system, which 
represents an interesting theoretical challenge for the future. Please note 
that a Gaussian approximation with 6^ = 2Tfc offers an alternative closure of 
moment equations for nonlinear friction function, but is unsuitable for the 
analysis of the disordered state at vanishing noise D as it would imply also 
a vanishing temperature, which for active particles is certainly not the case. 

Our kinetic equations can also be used to analyze the stability of the 
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spatially homogeneous state with respect to spatial perturbations. A corre- 
sponding instability in large scale collective motion has been predicted from 
the analysis of hydrodynamic equations of self-p r opelled particle sys tems 



fiBertin et al.l . I2006L l2009l : ISimha and Ramaswamvl . |2002b|ja|; llhld. 120111) and 



large density fluctuations were also reported in nu merical studies ( IGregoire and Chatd 



2OO4J : iGinelli et al.l . l2010t iPeruani et al.l . |2010[ ). Similar density inhomo 



geneities can be observed in numerical simulations of our model for e <^ L 
as shown in Fig. HI Thus, we should emphasize that our results for the 
homogeneous case hold strictly speaking only in the limit of global coupling 
and only as an approximation for local coupling, where the interaction range 
is not to small e in comparison to the system size L. 

The approach presented in this paper, can be applied also to other swarm- 
ing mode ls, such as, the "es c ape fc pursuit" model introduced recently by the 
authors ( iRomanczuk et al.l . 120091 ). The formulation of kinetic equations for 
interacting agents allow to address the question of how swarming affects pop- 
ulation dynamics, in particular the population dispersal, at ecological length 
scales. A Velocity-alignment interaction between individuals may result in 
long range ordered collective motion, which corresponds to large scale pop- 
ulation fluxes, and therefore the population dynamics can not be described 
by simple diffusion equations. 

We believe that the general approach of formulating individual based 
models in terms of stochastic differential equations and the derivation of 
macroscopic equations based on the corresponding nonlinear Fokker-Planck 
equation can be very useful to address many different problems of interacting 
agents in ecology. In particular, in the context of ecological modelling, it 
allows to establish a direct link between the individual based description 
and the macroscopic equation for systems of interacting agents. This link 
between the different levels of description is essential to gain a profound 
understanding of the self-organization on the population level but its direct 
mathematical derivation represents often a major challenge. 
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